#include "Mobius.h"
#include "Polynomicon.h"
#include "FourierTron.h"
#include <iostream>


using namespace std;
using namespace SeamStress;


Mobius::Mobius(Polynomicon* pol, vector<Seamstress*>* ss, unsigned long int nt, mpfr_prec_t p) : pins(this, ss), prec(p), reg(0,prec), reh(0,prec), retempval1(0,prec), retempval2(0,prec), retempval3(0,prec)
{
  if(nt==0){nthreads=1;}
  else{nthreads=nt;}
  if(nthreads > ss->size()){nthreads = ss->size();}
  
  poly = pol;
  
  g.set(0,0,prec);
  h.set(0,0,prec);
  tempval1.set(0,0,prec);
  tempval2.set(0,0,prec);
  tempval3.set(0,0,prec);
  
  threadvec1.assign(nthreads, tempval1);
  rethreadvec1.assign(nthreads, retempval1);
  
  tempvec1.assign(2, tempval1);
  retempvec1.assign(2, retempval1);
  
  h1.assign(2, tempvec1);
  h2.assign(2, tempvec1);
  htemp.assign(2, tempvec1);
  H.assign(2, tempvec1);
}


Mobius::~Mobius()
{
  
}


//takes inputpoly(x) of size N, transforms it to a polynomial in y, divided by (y+shift)^(N-1), for x = (a*y + b)/(c*y + d).  inputpoly gets sets to the coefficients of the numerator polynomial, and shift gets set to d/c
void Mobius::transform(const cxfloat& a, const cxfloat& b, const cxfloat& c, const cxfloat& d, vector<cxfloat>& inputpoly, cxfloat& shift)
{
  unsigned long int thread_count = nthreads;
  if(((inputpoly.size())>>1) < thread_count){thread_count = ((inputpoly.size())>>1);}
  if((inputpoly.size() - thread_count - 1)<thread_count){thread_count = (inputpoly.size() - thread_count - 1);}
  h.set(c);
  h.reciprocal();
  mul(a, h, g);
  mul(d, h, shift);
  h.sqr();
  mul(b, c, tempval1);
  mul(a, d, tempval2);
  tempval1.sub(tempval2);
  h.mul(tempval1);
  
  poly->shiftPoly(inputpoly, g);
  
  //now inputpoly is a poly in v = (x-g).  next we change variables to t = v/h
  poly->scalePoly(inputpoly, h);
  
  //now inputpoly is a poly in t.  next we change variables to u = 1/t, which just amounts to some swapping
  unsigned long int count = inputpoly.size() - 1;
  for(unsigned long int i=0;i<(((inputpoly.size())>>1)-1);i++)
  {
    __builtin_prefetch(&(inputpoly[i+1]), 1, 3);
    __builtin_prefetch(&(inputpoly[count-1]), 1, 3);
    swap(inputpoly[i], inputpoly[count]);
    count--;
  }
  swap(inputpoly[(((inputpoly.size())>>1)-1)], inputpoly[count]);
  
  //now we change variables from u to y, with y = u - shift
  poly->shiftPoly(inputpoly, shift);
}


void Mobius::transform(const refloat& a, const refloat& b, const refloat& c, const refloat& d, refloat_vec& inputpoly, refloat& shift)
{
  unsigned long int thread_count = nthreads;
  if(((inputpoly.size())>>1) < thread_count){thread_count = ((inputpoly.size())>>1);}
  if((inputpoly.size() - thread_count - 1)<thread_count){thread_count = (inputpoly.size() - thread_count - 1);}
  reh.set(c);
  reh.reciprocal();
  mul(a, reh, reg);
  mul(d, reh, shift);
  reh.sqr();
  mul(b, c, retempval1);
  mul(a, d, retempval2);
  retempval1.sub(retempval2);
  reh.mul(retempval1);
  
  poly->shiftPoly(inputpoly, reg);
  
  //now inputpoly is a poly in v = (x-g).  next we change variables to t = v/h
  poly->scalePoly(inputpoly, reh);
  
  //now inputpoly is a poly in t.  next we change variables to u = 1/t, which just amounts to some swapping
  unsigned long int count = inputpoly.size() - 1;
  for(unsigned long int i=0;i<(((inputpoly.size())>>1)-1);i++)
  {
    __builtin_prefetch(&(inputpoly[i+1]), 1, 3);
    __builtin_prefetch(&(inputpoly[count-1]), 1, 3);
    swap(inputpoly[i], inputpoly[count]);
    count--;
  }
  swap(inputpoly[(((inputpoly.size())>>1)-1)], inputpoly[count]);
  
  //now we change variables from u to y, with y = u - shift
  poly->shiftPoly(inputpoly, shift);
}


//calculate first nterms terms of power series of (x + l)^-M around x=0
void Mobius::poleSeries(const cxfloat& l, unsigned long int M, unsigned long int nterms, vector<cxfloat>& series)
{
  unsigned long int thread_count = nthreads;
  if(((nterms)>>1) < thread_count){thread_count = ((nterms)>>1);}
  if((nterms - thread_count - 1)<thread_count){thread_count = (nterms - thread_count - 1);}
  
  //set series[i] to (-1)^i * (x+l)^-(M+i)
  tempval1.set(l);tempval1.reciprocal();
  pow(M, tempval1, tempval2);
  if(series.size() < nterms)
  {
    for(unsigned long int i=0;i<series.size();i++){series[i].set(tempval2);}
    series.resize(nterms, tempval2);
  }
  else
  {
    series.resize(nterms, tempval2);
    for(unsigned long int i=0;i<series.size();i++){series[i].set(tempval2);}
  }
  tempval1.neg();
  poly->scalePoly(series, tempval1);
  
  //multiply series[i] by (M+i-1)!/((M-1)! * i!)
  poly->calculateFactorials(nterms + M);
  poly->factorialRecipVal(M-1, tempval1);
  thread_tot = nthreads;
  thread_var1 = M-1;
  tempvecpointer1 = &series;
  pins.sewStraight(&Mobius::poleSeriesThread1, thread_tot);
}


void Mobius::poleSeriesThread1(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  unsigned long int i = w;
  unsigned long int max = tempvecpointer1->size() - thread_tot;
  
  __builtin_prefetch(&((*tempvecpointer1)[i]), 1, 3);
  while(i < max)
  {
    __builtin_prefetch(&((*tempvecpointer1)[i+thread_tot]), 1, 3);
    poly->factorialVal(thread_var1 + i, threadvec1[w]);
    (*tempvecpointer1)[i].mul(threadvec1[w]);
    (*tempvecpointer1)[i].mul(tempval1);
    poly->factorialRecipVal(i, threadvec1[w]);
    (*tempvecpointer1)[i].mul(threadvec1[w]);
    
    i += thread_tot;
  }
  poly->factorialVal(thread_var1 + i, threadvec1[w]);
  (*tempvecpointer1)[i].mul(threadvec1[w]);
  (*tempvecpointer1)[i].mul(tempval1);
  poly->factorialRecipVal(i, threadvec1[w]);
  (*tempvecpointer1)[i].mul(threadvec1[w]);
}


void Mobius::poleSeries(const refloat& l, unsigned long int M, unsigned long int nterms, refloat_vec& series)
{
  unsigned long int thread_count = nthreads;
  if(((nterms)>>1) < thread_count){thread_count = ((nterms)>>1);}
  if((nterms - thread_count - 1)<thread_count){thread_count = (nterms - thread_count - 1);}
  
  //set series[i] to (-1)^i * (x+l)^-(M+i)
  retempval1.set(l);retempval1.reciprocal();
  pow(M, retempval1, retempval2);
  if(series.size() < nterms)
  {
    for(unsigned long int i=0;i<series.size();i++){series[i].set(retempval2);}
    series.resize(nterms, retempval2);
  }
  else
  {
    series.resize(nterms, retempval2);
    for(unsigned long int i=0;i<series.size();i++){series[i].set(retempval2);}
  }
  retempval1.neg();
  poly->scalePoly(series, retempval1);
  
  //multiply series[i] by (M+i-1)!/((M-1)! * i!)
  poly->calculateFactorials(nterms + M);
  poly->factorialRecipVal(M-1, retempval1);
  thread_tot = thread_count;
  thread_var1 = M-1;
  retempvecpointer1 = &series;
  pins.sewStraight(&Mobius::poleSeriesThreadReal1, thread_tot);
}


void Mobius::poleSeriesThreadReal1(void* arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  unsigned long int i = w;
  unsigned long int max = retempvecpointer1->size() - thread_tot;
  
  __builtin_prefetch(&((*retempvecpointer1)[i]), 1, 3);
  while(i < max)
  {
    __builtin_prefetch(&((*retempvecpointer1)[i+thread_tot]), 1, 3);
    poly->factorialVal(thread_var1 + i, rethreadvec1[w]);
    (*retempvecpointer1)[i].mul(rethreadvec1[w]);
    (*retempvecpointer1)[i].mul(retempval1);
    poly->factorialRecipVal(i, rethreadvec1[w]);
    (*retempvecpointer1)[i].mul(rethreadvec1[w]);
    
    i += thread_tot;
  }
  poly->factorialVal(thread_var1 + i, rethreadvec1[w]);
  (*retempvecpointer1)[i].mul(rethreadvec1[w]);
  (*retempvecpointer1)[i].mul(retempval1);
  poly->factorialRecipVal(i, rethreadvec1[w]);
  (*retempvecpointer1)[i].mul(rethreadvec1[w]);
}


void Mobius::transformDers(const cxfloat& a, const cxfloat& b, const cxfloat& c, const cxfloat& d, vector<cxfloat>& inputpoly, const cxfloat& shift)
{
  transform(a,b,c,d, inputpoly, tempval3);
  
  poly->shiftPoly(inputpoly, shift);
  tempval3.add(shift);
  
  poleSeries(tempval3, inputpoly.size() - 1, inputpoly.size(), tempvec1);
  poly->multiplyPoly(inputpoly, tempvec1, tempvec2);
  tempvec2.resize(inputpoly.size(), tempval1);
  
  unsigned long int M = inputpoly.size()-1;
  for(unsigned long int i=0;i<(M);i++)
  {
    __builtin_prefetch(&(inputpoly[i+1]), 1, 3);
    __builtin_prefetch(&(tempvec2[i+1]), 0, 3);
    inputpoly[i].set(tempvec2[i]);
  }
  inputpoly[M].set(tempvec2[M]);
}


void Mobius::transformDers(const refloat& a, const refloat& b, const refloat& c, const refloat& d, refloat_vec& inputpoly, const refloat& shift)
{
  transform(a,b,c,d, inputpoly, retempval3);
  
  
  poly->shiftPoly(inputpoly, shift);
  retempval3.add(shift);
  
  poleSeries(retempval3, inputpoly.size() - 1, inputpoly.size(), retempvec1);
  poly->multiplyPoly(inputpoly, retempvec1, retempvec2);
  retempvec2.resize(inputpoly.size(), retempval1);
  
  unsigned long int M = inputpoly.size()-1;
  for(unsigned long int i=0;i<(M);i++)
  {
    __builtin_prefetch(&(inputpoly[i+1]), 1, 3);
    __builtin_prefetch(&(retempvec2[i+1]), 0, 3);
    inputpoly[i].set(retempvec2[i]);
  }
  inputpoly[M].set(retempvec2[M]);
}


void Mobius::matrixMul2D(const vector<vector<cxfloat> >& M1, const vector<vector<cxfloat> >& M2, vector<vector<cxfloat> >& M3)
{
  for(unsigned long int i=0;i<2;i++)
  {
    for(unsigned long int j=0;j<2;j++)
    {
      mul(M1[i][0], M2[0][j], tempval1);
      mul(M1[i][1], M2[1][j], M3[i][j]);
      M3[i][j].add(tempval1);
    }
  }
}


void Mobius::matrixInv2D(const vector<vector<cxfloat> >& M1, vector<vector<cxfloat> >& M2)
{
  mul(M1[0][0], M1[1][1], tempval1);
  mul(M1[0][1], M1[1][0], M2[0][0]);
  tempval1.sub(M2[0][0]);
  tempval1.reciprocal();
  
  mul(M1[1][1], tempval1, M2[0][0]);
  mul(M1[0][0], tempval1, M2[1][1]);
  mul(M1[0][1], tempval1, M2[0][1]);M2[0][1].neg();
  mul(M1[1][0], tempval1, M2[1][0]);M2[1][0].neg();
}


void Mobius::map(const cxfloat& z1, const cxfloat& z2, const cxfloat& z3, const cxfloat& w1, const cxfloat& w2, const cxfloat& w3)
{
  sub(z2, z3, h1[0][0]);
  sub(z2, z1, h1[1][0]);
  mul(z1, h1[0][0], h1[0][1]);h1[0][1].neg();
  mul(z3, h1[1][0], h1[1][1]);h1[1][1].neg();
  
  sub(w2, w3, h2[0][0]);
  sub(w2, w1, h2[1][0]);
  mul(w1, h2[0][0], h2[0][1]);h2[0][1].neg();
  mul(w3, h2[1][0], h2[1][1]);h2[1][1].neg();
  
  
  matrixInv2D(h2, htemp);
  matrixMul2D(htemp, h1, H);
  
  matrixInv2D(H, h2);
}


void Mobius::compose(const cxfloat& a, const cxfloat& b, const cxfloat& c, const cxfloat& d, const cxfloat& A, const cxfloat& B, const cxfloat& C, const cxfloat& D, cxfloat& a2, cxfloat& b2, cxfloat& c2, cxfloat& d2, cxfloat& A2, cxfloat& B2, cxfloat& C2, cxfloat& D2)
{
  htemp[0][0].set(a);
  htemp[0][1].set(b);
  htemp[1][0].set(c);
  htemp[1][1].set(d);
  
  h1[0][0].set(a2);
  h1[0][1].set(b2);
  h1[1][0].set(c2);
  h1[1][1].set(d2);
  
  matrixMul2D(htemp, h1, h2);
  a2.set(h2[0][0]);
  b2.set(h2[0][1]);
  c2.set(h2[1][0]);
  d2.set(h2[1][1]);
  
  
  htemp[0][0].set(A);
  htemp[0][1].set(B);
  htemp[1][0].set(C);
  htemp[1][1].set(D);
  
  h1[0][0].set(A2);
  h1[0][1].set(B2);
  h1[1][0].set(C2);
  h1[1][1].set(D2);
  
  matrixMul2D(h1, htemp, h2);
  A2.set(h2[0][0]);
  B2.set(h2[0][1]);
  C2.set(h2[1][0]);
  D2.set(h2[1][1]);
}




